#Replication file for "The Logic of Ambitious Legislators in Fluid Party Systems
#Legislative Studies Quarterly
#Thomas Mustillo, John Polga-Hecimovich, Francisco Sanchez
#June 2024

# Packages
library("jtools")
library("gridExtra")
library("tidyverse")
library("viridis")
library("patchwork")
library("hrbrthemes")
library("circlize")
library("stargazer")
library("distributional")
library("ggdist")
library("ggthemes")
library("broom.mixed")
library("broom")
library("ggplot2")
library("dplyr")
library("rms")
library("marginaleffects")
library("modelsummary")
library("readxl")

###############
############### Descriptive graphics
###############

###############
############### Graphic 1
###############
# Load dataset
sankey <- read.csv (file = 'Sankey.csv', row.names = 1)
View(sankey)
library(networkD3)

# Long format
data_long <- sankey %>%
  rownames_to_column %>%
  gather(key = 'key', value = 'value', -rowname) %>%
  filter(value > 0)
colnames(data_long) <- c("source", "target", "value")
data_long$target <- paste(data_long$target, " ", sep="")

# Create a node data frame: it lists every entities involved in the flow
nodes <- data.frame(name=c(as.character(data_long$source), as.character(data_long$target)) %>% unique())

# With networkD3, connection must be provided using id, not using real name like in the links dataframe. So we need to reformat it.
data_long$IDsource=match(data_long$source, nodes$name)-1 
data_long$IDtarget=match(data_long$target, nodes$name)-1

# prepare color scale
# ColourScal ='d3.scaleOrdinal() .range(["#FDE725FF","#B4DE2CFF","#6DCD59FF","#35B779FF","#1F9E89FF","#26828EFF","#31688EFF","#3E4A89FF","#482878FF","#440154FF"])'
ColourScal ='d3.scaleOrdinal(d3.schemeCategory20)'

# Make the Network. Graphic 1
p <- sankeyNetwork(Links = data_long, Nodes = nodes,
                   Source = "IDsource", Target = "IDtarget",
                   Value = "value", NodeID = "name", 
                   sinksRight=FALSE, colourScale=ColourScal, nodeWidth=20, fontSize=10, nodePadding=6)

## save to html
library(htmlwidgets)
htmlwidgets::saveWidget(widget = p, file = "Sankey1.html")
# save html to png
webshot::webshot(url = "Sankey1.html",
                 file = "Sankey1.png" )
png('Sankey1.png', width = 400,height = 600)
print(p)
dev.off()


###############
############### Graphic 2
############### 
sankey2 <- read.csv (file = 'Sankey2.csv', row.names = 1)
View(sankey2)
library(networkD3)

# Long format
data_long2 <- sankey2 %>%
  rownames_to_column %>%
  gather(key = 'key', value = 'value', -rowname) %>%
  filter(value > 0)
colnames(data_long2) <- c("source", "target", "value")
data_long2$target <- paste(data_long2$target, " ", sep="")

# Create a node data frame: it lists every entity involved in the flow
nodes2 <- data.frame(name=c(as.character(data_long2$source), as.character(data_long2$target)) %>% unique())

# With networkD3, connection must be provided using id, not using real name like in the links dataframe. So we need to reformat it.
data_long2$IDsource=match(data_long2$source, nodes2$name)-1 
data_long2$IDtarget=match(data_long2$target, nodes2$name)-1

# prepare color scale
ColourScal ='d3.scaleOrdinal(d3.schemeCategory20)'

# Make the Network
sankeyNetwork(Links = data_long2, Nodes = nodes2,
              Source = "IDsource", Target = "IDtarget",
              Value = "value", NodeID = "name", 
              sinksRight=FALSE, colourScale=ColourScal, nodeWidth=20, fontSize=10, nodePadding=6)

## save to html
library(htmlwidgets)
htmlwidgets::saveWidget(widget = q, file = "Sankey2.html")
# save html to png
webshot::webshot(url = "Sankey2.html",
                 file = "Sankey2.png" )
png('Sankey2.png', width = 400,height = 600)
print(p)
dev.off()



###############
############### Public records data, tables, and figures
###############

#Load data
switching <- read_excel("Ecuador_PublicData.xlsx")

#Model 1: All years
    model1 <- glm(switch ~ 
      v2pariglef +
      v2xpa_popul +
      oldparty_ideology +
      oldparty_popul +
      pres_party*approval_smoothed +
      pres_top2 +
      dm + 
      total_terms +
      alliance +
      woman,
      data=switching,
      family = binomial(link = "logit"))
#coefficient map
cm <- c("v2pariglef" = "Party ideology", "v2xpa_popul" = "Party populism", "oldparty_ideology" = "Previous party ideology", "oldparty_popul" = "Previous party populism", "pres_party" = "President's party", "approval_smoothed" = "Presidential approval", "pres_party:approval_smoothed" = "Pres's party * Approval", "pres_top2" = "Presidential coattails", "dm" = "District magnitude", "total_terms" = "Total terms", "alliance" = "Alliance", "woman" = "Female", "(Intercept)" = "Intercept")
    #table reporting odds ratios (rather than log odds ratios)
    modelsummary(model1, exponentiate = TRUE, coef_map = cm)
    
#Model 5: Subset 1, pre-1996
model5 <- glm(switch ~ 
      v2pariglef +
      v2xpa_popul + 
      oldparty_ideology + 
      oldparty_popul + 
      pres_party*approval_smoothed + 
      pres_top2 + 
      dm  + 
      total_terms + 
      alliance + 
      woman,
      data=switching,
      family = binomial(link = "logit"),
      subset = year < 1996)
    cm <- c("v2pariglef" = "Party Ideology", "v2xpa_popul" = "Party Populism", "oldparty_ideology" = "Previous Party Ideology", "oldparty_popul" = "Previous party populism", "pres_party" = "President's Party", "approval_smoothed" = "Presidential Approval", "pres_party:approval_smoothed" = "PresParty*Approval", "pres_top2" = "Presidential Coattails", "dm" = "District Magnitude", "total_terms" = "Total Terms", "alliance" = "Alliance", "woman" = "Woman")
    modelsummary(model5, exponentiate = TRUE, coef_map = cm)
    
#Model 7: Subset 2, post-1996
    model7 <- glm(switch ~ 
      v2pariglef  + 
      v2xpa_popul + 
      oldparty_ideology + 
      oldparty_popul + 
      pres_party*approval_smoothed + 
      pres_top2 + 
      dm + 
      total_terms + 
      alliance + 
      woman,
      data=switching,
      family = binomial(link = "logit"),
      subset = year >= 1996)
    cm <- c("v2pariglef" = "Party Ideology", "v2xpa_popul" = "Party Populism", "oldparty_ideology" = "Previous Party Ideology", "oldparty_popul" = "Previous party populism", "pres_party" = "President's Party", "approval_smoothed" = "Presidential Approval", "pres_party:approval_smoothed" = "PresParty*Approval", "pres_top2" = "Presidential Coattails", "dm" = "District Magnitude", "total_terms" = "Total Terms", "alliance" = "Alliance", "woman" = "Woman")
    modelsummary(model7, exponentiate = TRUE, coef_map = cm)
    
    models <- list(
      "Full sample" = model1,
      "Pre-1996" = model5,
      "Post-1996" = model7
    )
    cm <- c("v2pariglef" = "Party ideology", "v2xpa_popul" = "Party populism", "oldparty_ideology" = "Previous party ideology", "oldparty_popul" = "Previous party populism", "pres_party" = "President's party", "approval_smoothed" = "Presidential approval", "pres_party:approval_smoothed" = "Pres's party * Approval", "pres_top2" = "Presidential coattails", "dm" = "District magnitude", "total_terms" = "Total terms", "alliance" = "Alliance", "woman" = "Female", "(Intercept)" = "Intercept")
    #to view
    modelsummary(models, exponentiate = TRUE, stars = TRUE, coef_map = cm, notes = list('Coefficient estimates are odds ratios.'))
    #to tex
    modelsummary(models, output = "switching.tex", exponentiate = TRUE, stars = TRUE, coef_map = cm, notes = list('Coefficient estimates are odds ratios.'))
    
    pnew_ideo <- plot_predictions(model1, condition = "v2pariglef", rug = TRUE) + 
      labs(x = "Ideology of current party", y = "Pr(Switching)") +
      theme_bw() +
      expand_limits(x = c(-3, 3)) +
      ggtitle("Panel C")
    
    pnew_pop <- plot_predictions(model1, condition = "v2xpa_popul", rug = TRUE) + 
      labs(x = "Populism of current party", y = "Pr(Switching)") +
      theme_bw() +
      expand_limits(x = c(0, 1)) +
      ggtitle("Panel D")
    
    pold_ideo <- plot_predictions(model1, condition = "oldparty_ideology", rug = TRUE) + 
      labs(x = "Ideology of former party", y = "Pr(Switching)") +
      theme_bw() +
      expand_limits(x = c(-3, 3)) +
      ggtitle("Panel A")
    
    pold_pop <- plot_predictions(model1, condition = "oldparty_popul", rug = TRUE) + 
      labs(x = "Populism of former party", y = "Pr(Switching)") +
      theme_bw() +
      expand_limits(x = c(0, 1)) +
      ggtitle("Panel B")
    
    grid.arrange(pold_ideo, pold_pop, pnew_ideo, pnew_pop, nrow = 2)
    plot <- grid.arrange(pold_ideo, pold_pop, pnew_ideo, pnew_pop, nrow = 2)
    
  
###############
############### PELA data and tables
###############

#Load data
pela <- read_excel("Ecuador_PELA.xlsx")
pela <- as.data.frame(pela)

#Data management
  # Create dummy variables for the factor in "dis3"
  # En ocasiones sucede que un Diputado que ha sido elegido en las listas de un determinado partido 
    #decide posteriormente abandonar su vinculación con ese partido político. 
    #¿qué cree Ud. que debería hacer el Diputado
    Retain <- ifelse(pela$dis3 == '1', 1, 0) # Conservar su escaño e integrarse en otro Grupo Parlamentario

    # Create dummy variables for the factors in "rep5"
    Department <- ifelse(pela$rep5 == '2', 1, 0) # A todos los electores de su Departamento

    pela <- cbind(pela, Retain, Department)

  #Model 1: 1994-2017 (all six waves)
    pela1 <- glm(switch ~
      id1 +
      id2 + 
      Retain +
      rep1,
      data=pela,
      family = binomial(link = "logit"))
    modelsummary(pela1, exponentiate = TRUE)

  #Model 2: 1994-2017 (all six waves)
    pela2 <- glm(switch ~
      id1 +
      id2 + 
      dis2 +
      Retain +
      rep1 +
      Department, 
      data=pela,
      family = binomial(link = "logit"))
    modelsummary(pela2, exponentiate = TRUE)

  #Model 3: 1994-2017 (all six waves)
    pela3 <- glm(switch ~
      id1 +
      id2 + 
      dis2 +
      Retain +
      rep1 +
      Department +
      inst102 +
      inst107,
      data=pela,
      family = binomial(link = "logit"))
    cm <- c("id1" = "Ideology", "id2" = "Party ideology", "dis2" = "Expel dep. against discipline", "Retain" = "Retain", "rep1" = "Resources province", "Department" = "Provincial rep.", "inst102" = "Confidence parties", "inst107" = "Confidence personal activ.", "(Intercept)" = "Intercept")
    modelsummary(pela3, exponentiate = FALSE, coef_map = cm)

  #All models combined
    models <- list(
      "(1)" = pela1,
     "(2)" = pela2,
     "(3)" = pela3
    )

    cm <- c("id1" = "Self Ideology", "id2" = "Party ideology", "dis2" = "Expel dep. against discipline", "Retain" = "Retain Seat", "rep1" = "Resources province", "Department" = "Provincial rep.", "inst102" = "Confidence parties", "inst107" = "Confidence personal activ.", "(Intercept)" = "Intercept")
    #to view
    modelsummary(models, exponentiate = TRUE, stars = TRUE, coef_map = cm, notes = list('Coefficient estimates are odds ratios.'))
    #to tex
    modelsummary(models, output = "pela.tex", exponentiate = TRUE, stars = TRUE, coef_map = cm, notes = list('Coefficient estimates are odds ratios.'))


